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ABSTRACT 

In this paper, we investigate the structures of galaxies which either have or 
have had three BHs using iV-body simulations, and compare them with those of 
galaxies with binary BHs. We found that the cusp region of a galaxy which have 
(or had) triple BHs is significantly larger and less dense than that of a galaxy 
with binary BHs of the same mass. Moreover, the size of the cusp region depends 
strongly on the evolution history of triple BHs, while in the case of binary BHs, 
the size of the cusp is determined by the mass of the BHs. In galaxies which have 
(or had) three BHs, there is a region with significant radial velocity anisotropy, 
while such a region is not observed in galaxies with binary BH. These differences 
come from the fact that with triple BHs the energy deposit to the central region of 
the galaxy can be much larger due to multiple binary-single BH scatterings. Our 
result suggests that we can discriminate between galaxies which experienced triple 
BH interactions with those which did not, through the observable signatures such 
as the cusp size and velocity anisotropy. 

Subject headings: black hole: physics — black hole: binary — galaxies: nuclei 
— galaxies: structure — stellar dynamics — Three-body problem:general — 
methods: n-body simulations 



Introduction 



According to recent observations, elliptical galaxies can be classified into two g roups: 
"weak-cusp" galaxies and "strong-cusp" galaxies (ILauer et al.lll995l : iFaber et al.lll997l ). The 
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central surface brightness profiles of the weak-cusp galaxies are expressed as S(-R) oc i? -7 
with 7 < 0.3, and those of the strong-cusp galaxy the same formula with 7 > 0.5. The slope 
of the volume density profile of strong-cusp galaxies is around —2, being consistent with the 
isothermal cusp. Such a steep cusp is naturally formed in dissipative process involving gas 
dynamics and star formation. On the other hand, the weak cusp corresponds to the slope of 
the volume density shallower than —1, which is not likely to be formed through dissipative 
process. 

One possible way to form a weak cusp is the merging of two galaxies containing black 

holes (BHs). When two galaxies, each with a black hole, merge, two BHs s ink toward the cen- 

terof the merger remnant to form a binary through the dynamical friction (IBegelman. Blandford. fc Rees 
19801 ). The back reaction of the dynamical f riction heats up field stars. As a result, a shallow 



cusp of stars develops in the central region (lEbisuzaki et al.lll99ll ; iNakano fc Makino 111999 
MerrittlliooiSl ). 



One problem with this binary B H sce nario is what would be the final fate of the binary 
BH. IBegelman. Blandford. fc Reed (jl980j) pointed out that the merging timescale of the 
binary BH might be very long, after the binary BH ejected out the stars which can interact 
with the binary (loss cone depletion). The stars will be supplied in the timescale of the 
relaxation time, which is much long er than the Hubble tim e . Recent iV-body simulations 
confi rmed this theoretical estimate (IMakino fc Funatol 12004 ; I Berczik. Merritt. fc Spurzem 
2005h . 



If galaxies are formed through hierarchical clusterings, in many cases, a binary BH is 
formed after a merger event. If there were sufficient gas left in the merger, interaction with 
gas might lead to the quick merging of the two BHs. However, in the case of "dry" mergers 
which would result in the formation of giant ellipticals, by definition not much gas is left 
and it would be difficult for two BHs to merge. 

If one galaxy with binary BH and the other with single BH merge, the central BHs form a 
triple system. Iwasawa, Funato & Makino (2006, hereafter referred to as Paper I) investigated 
the evolution of triple BH system in the galactic center, using iV-body simu lations. They 
found that the strong binary-single BH interact ion (IMakino fc Ebisuzakill 19941 ) and the Kozai 
cycle ( jKozaill 19621 ; iBlas. Lee, fc Socrates Il2002l ) drives the eccentricity of t he BH binary high 



enoug h that two BHs merge quickly through gravitational wave radiation. iHoffman fc Loeb 



( 120071 ) performed statistical simulations of evolution of central BH systems and reached a 
similar conclusion. 



This paper is a follow-up of Paper I. In this paper, we investigate the structure of 
a galaxy containing three BHs. We also investigated their observational properties which 



- 3- 



would help us to find galaxies which have (or had) triple BHs. 

We performed iV-body simulations of the evolution of triple (or binary) BH in a host 
galaxy, in order to study the dynamical evolution of the structure of stellar systems contain- 
ing BHs. In our simulations, both dynamical evolution of BHs and that of field stars are 
integrated consistently. The interaction between BHs affects not only the evolution of them- 
selves but also the spatial and kinematic structure of field stars around them. In turn, the 
distribution of field stars affects the interaction between BHs. To understand the structure 
of galaxies containing BHs, a self-consistent simulation in which the orbits of BHs and stars 
are treated self-consistently is essential. 

The structure of this paper is as follows. In section 2, we describe the initial models and 
the method of our numerical simulations. In section 3, we show the effect of the BH triples 
dynamics on the structure of the galaxy. Summary and discussion are given in section 4. 



Initial Models and Numerical Methods 



2.1. Initial Models 



The initial setup of the models is basically the same as that in Paper I. For the galaxy 
model, we used a King model with Wn = 11, where Wo is th e nondimensional central po- 
tential of King models (Kind 1966; iBinnev &: Tremaine 1987 ). We adopted the standard 
iV-body units (IHeggie fc Mathieulll986l ). in which m ga i = 1.0,E gal = — 0.25, G = 1. Here, 
m g ai and E ga i are the total mass and total binding energy of the galaxy not including BHs 
and G is the gravitational constant. We placed three equal- mass BH particles in iV-body 
models of a spherical galaxy. Two of three BHs particles are initially placed at the po- 
sition (±0.005,0.0,0.0) with velocity (0.0, ±0.15, 0.0), and the third one is at the position 
(0.1, 0.0, 0.0) with velocity (0.0, 0.0, 0.0). In order to compare the structure of a galaxy with 
three BHs with that of a galaxy with two BHs, we also performed simulations of galaxies 
with two BHs. Total mass of BHs in two-BH models was kept to be the same as that in 
the corresponding three-BH models. For two-BH model, two BHs are placed at the posi- 
tion (±0.1, 0.0, 0.0) with velocity (0.0, ±0.1, 0.0). The quantitative properties of models and 
initial conditions of our simulation is summarized in table [TJ 

We set the number of FS particles N-ps — 262144 in all simulations. In order to inves- 
tigate the effect of triple BHs on the galaxy, we set the mass ratio between galaxy and total 
mass of BHs to, mBH.tot/^gai ~ 0.003, where m B H,tot is total mass of BH, in our standard 



model. This value is what is suggested by rece nt observations (IKormendy fc Richstondll995 



Magorrian et al.lll998l ; iMarconi fc Huntll2003l ). Table [2] shows all models 
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2.2. Numerical Method 

The numerical method is the same as that we used in Paper I. To summarize, we used 
a softened Newton gravity for the forces between field stars and those between each BH 
particle and field stars, while for forces between BH particles we applied post-Newtonian 
approximation to include the back reaction of gravitational radiation to the BHs. For the 
term corresponding to the radiation o f the gravitati onal wave, we used approximately 2.5- 



order post-Newtonian approximation (iDamourl 119871 ) . 

We used the Plummer softening and set the value eps-FS — 10~ 4 , cfs-bh = 10~ 7 and 
cbh-bh = 0.0, where exx-YY is softening parameters used for a XX part icle - a YY particle 



inter action. Time integration scheme is the 4-th order Hermite scheme ( IMakino fe Aarseth 



19921 ) with individual variable time steps. 



In order to calculate the acceleration due to field particles, we used GRAPE-6 ( IMakino et al 



20031 ). the special-purpose computer for the gravitational A^-body problem. Forces from BH 
particles, both Newtonian and gravitational wave terms, are calculated on the host com- 
puter. In all runs, the energy (corrected for the loss of energy through GW radiation) is 
conserved better than 0.2 % 



2.3. Physical Scales 

The mass unit in our simulation corresponds to 10 n M & . Thus, the mass of each BH 
particle is 1O 8 M in model MIT. We assume that the velocity dispersion of the initial King 
model at the King radius corresponds to 300 km/s. In other words, we set the light velocity 
c in the A^-body unit to 1006. The unit of length and time are about 4.9kpc and 16Myr, 
respectively. The king radius is about lOOpc. 



3. Result 

3.1. Evolution of BH System 

Figure [1] shows the evolution of semi-major axis a and eccentricity e of BH binaries in 
all models. For models with triple BHs (model name MxT), the definition of the binary 
of BHs is the most strongly bound pair of three BHs. For the model MIT, jumps in a 
and e at T ~ 0.45, 0.96, 2.2, 3.7 and 4.8 are results of interactions between the BH binary 
and the third BH. At T ~ 4.8, two BHs merged. This merging is driven by the impulsive 
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Table 1. Model Parameters 



Parameter 


Symbol 


Value 


Mass of galaxy 


m ga i 


l.U 


Mass of FS 


raps 


1/262144 


Mass of BH 




0.001 -0.0045 


Total mass of BHs 


^BH,tot 


0.003 - 0.009 


Number of FS 


N FS 


262144 


Number of BH 


N B n 


2or3 


Gravitational constant 


G 


1 


Total energy 


E ga \ 


-0.25 


Softening between BHs 


CBH-BH 


0.0 


Softening between field stars and BH 


^fs-bh 


io- 7 


Softening between field stars 


cfs-fs 


0.001 



Table 2. Model List 



Model name 




A^BH 


Wo 


^BR- 


final state 


MIT 


262144 


3 


11 


0.001 


merge 


M2T 


262144 


3 


11 


0.002 


merge 


M3T 


262144 


3 


11 


0.003 


merge 


M1B 


262144 


2 


11 


0.0015 


no merge 


M2B 


262144 


2 


11 


0.003 


no merge 


M3B 


262144 


2 


11 


0.0045 


no merge 
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eccentricity growth caused by a binary-single interaction (Paper I), followed by the spiral-in 
due to gravitational wave radiation. The merging criterion we used here is that the distance 
between two BHs becomes smaller than three times of the Schwarzschild radius of a BH. We 
replaced the merged BHs by a new single BH. This new BH has the same mass, position 
and velocity as those of the barycenter of original binary. 

As shown in Figure [TJ the merged BH and the third one form a binary in model MIT. 
In models M1B and MIT, after T = 4.8, a and e change slowly, and any BH merger does 
not occur until T = 10. In both models, a binary BH is left around the center of the host 
galaxy. Note that this slow shrink of the BH binary orbit is due to the refilling of the loss 
cone through two-body relaxation, which is enhanced because of the relatively small number 
of particles we used (iVps = 262144). In real galaxies the timescale of orbital evolution of 
binary BH would be much longer after the loss cone is depleted. 

In model M2T, the overall behavior of BHs is rather similar to that in model MIT. After 
strong triple interactions at T = 0.35 and 2.2, the binary merged through gravitational wave 
radiation. On the other hand, in model M3T, after the first triple interaction at T = 0.3, 
the binary BH hardens through the interaction with field stars and finally merged through 
gravitational wave radiation. The difference between these three models is at least partly 
just a coincidence, since the final result depends on the outcome of the first triple encounter 
at time around 0.4. 

On the other hand, the behavior of models with binary BHs are all very similar, showing 
slow hardening driven by the loss-cone refill through two-body relaxation. 



3.2. Density Structure 

3.2.1. Density Profile 

In Figure [2], we show the spatial (left) and surface (right) density profiles of models M3T 
and M3B at T = 3 and T = 6. In both models the central density decreases and a weak 
cusp develops in the central region by T = 3. They do not change significantly after T = 3. 
This cusp is expressed as p oc r~ 7 , where 7 ~ 0.5. For both m odels the slope of cusps can 



be explained by the simple theory of iNakano fc Makino I (Il999l ). 



From figure [2] we can see that there is a clear difference between models MIT and M1B. 
The size of the cusp region is bigger for MIT and the density is lower. This result indicates 
that a triple BH system is more efficient in heating the central region and ejecting the stars 
than a BH binary is. 
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In the case of a galaxy with a binary BH, a BH binary hardens in a monotonic fashion 
through interactions with FSs. In a galaxy with a triple BH system, a BH binary experiences 
multiple interactions with a single BH. In each event a single BH and the center of mass 
motion of binary BH acquires energy. This energy is then transferred to field stars through 
dynamical friction. As a result, the field stars are heated up and the cusp becomes wider. 
This mechanism is shown in Figure [3] more clearly. 



3.2.2. Evolution of Cusp Region 

In figures [3] and HI time evolution of the Lagrangian radii of field stars, the cusp radius 
and radial positions of BHs for all models are plotted. 



Following ICasertano fc Hutl (119851 ) . we defined the cusp radius as 

= piTi (i) 

1 cusp — ' \ J 

l^ipi 

where pi is the local mass density around the field particle with index % and the summation 
is done over field particles. The local mass density is defined as 

Pi = 4 3~> ( 2 ) 

where, m^i is the total FS particle mass contained within tq^ which is the distance from the 
field particle i to its sixth nearest neighbor. 

In top panel of figure [31 the evolution of model MIT is shown. We can see that the 
Lagrangian radii and cusp radius change rather impulsively at the same time as the strong 
triple scattering events occur and BHs are ejected out of the central region. This expansion 
is because of the indirect heating due to the removal of massive BH particles from the center. 
After the scattering events, BHs sink towards the center of the system through dynamical 
friction, and the Lagrangian radii show small expansion. Thus, the expansion is driven by 
the BH particles. In the case of model M1B, the expansion is smooth and smaller than that 
in model MIT. 

Figure 0] shows the evolution of the cusp radius and the mean cusp density of host 
galaxies for all models. We can clearly see that the expansion of the cusp radius in models 
MxT are driven by triple scattering events. In the case of model M3T, there was only one 
scattering event and therefore only one rapid expansion event for the cusp radius. In this 
model, the third BH is ejected out of the galaxy and never returned. Thus, one merged BH 
particle is left at the center of the galaxy. Slow decrease of the cusp radius is due to the 
thermal evolution of the field star system. 
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Figure [5] shows the density profile for all runs at T = 6. The density profile for all 
models are similar in their shapes. The density profile of model M3T is a bit steeper than 
those of other models. This is because only one binary-single scattering event is not enough 
to make the density profile p oc r~°- 5 . 

For models with triple BHs (MxT), the radius and density of the cusp are determined 
not only by the total mass of BHs, but also by the number of binary-single BH scattering 
events. On the other hand, in the case of binary BHs, the total mass of BHs determines 
the structure of the cusp. In addition, in all cases the cusps in models MxT are bigger than 
those in models MxB of the corresponding BH mass. 



3.3. Mass deficit 



The discussion in the previous subsection can be re-interpreted in terms of the so-called 
"mass deficit" (Milosavljevic & Merritt 2001, hereafter referred to as MM2001 ). They argued 
that one can estimate the mass removed by BH binary activity by comparing the singular 
isothermal profile and an actual cusp profile, and they argued that the difference, w hich they 
named "mass deficit" , retains the memory of the merging history of the black hole. I Graham 



( 120041 ) argued that MM2001's method to estimate the mass deficit would give too large 
values like md e f ~ lOm-BH , where is the mass deficit, while a more reasonable method 



would give ~ m-Bn- iMerrittl (120061 ) claimed that in repeated mergers m^ei ~ A'wibh, 



where N is the number of merging events, though his simulation result seems to suggest 
that such linear scaling does not hold for N > 3. In addition, simulations of the repeated 
merger of IMerrittl (I200q) ignored the fac t that additional BHs in reality comes with their 
host galaxies. iMakino fc Ebisuzakil (119961 ) have already shown that in the case of hierarchical 
repeated mergings, the density profile of merger remnant converge to a universal profile. In 
other words, mdef/mBH does not tell much about the past merger history. 

Figure M shows the ratio of mass deficit to total mass of BHs for all models. Here, we 
define the mass deficit as difference of mass at any given time with original mass within 
"break radius", defined as the radius at which the initial profile and that at the current 
time intersect with each other (see Figure [2]). In the case of binary BH runs (model MxB), 
^def/ m BH,tot reaches about unity by T ~ 1, and the increase after that is due to loss-cone 
refilling by relaxation. Therefore, in the case where relaxation is negligible, mj e f ~ ?^BH,tot- 
In the case of triple BH system, however, m^ef /^BH.tot does tell about the past history of 
the evolution of the triple system. Thus, a galaxy with unusually large rridef /iriBK,tot is a 
good candidate for the host of triple BH. 
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In Figure 7, we compare the two methods to estimate the mass deficit, one of which 
is our definition and the other is used by MM2001. The estimate mass from MM2001 
is systematically higher than real one. This difference is simply due to the fact that the 
method of MM2001 overestimates the initial mass with the break radius. We tried to use 



the method proposed by iGrahaml (12004 ). but it turned out that the fit by a Sersic profile is 



not appropriate for our galaxy model. 



3.4. Velocity Structure 

Figure [S] shows the velocity dispersion (top) and line-of-sight velocity (bottom) profiles 
at T = 6. In the central region, the velocity dispersion is increasing toward the center. Since 
the potential of BHs dominates the potential in the central region, the orbits of stars are 
almost Keplerian. The line-of-sight velocity dispersion shows the similar rise at the center. 

In Figure M, the velocity anisotropy for models MIT and M1B are compared. Here the 
anisotropy parameter is defined as 



where (v%) and (v%) are the tangential and radial mean square velocities. For model M1B, 
orbits are circular in the central region and f3 increases towards as r increases at T = 3 
and 6. From figures [3] and [9j it is clear that, after two BHs settled and formed a binary in 
the center, the orbits of stars in the cusp becomes tangentially anisotropic and in the outer 
region the orbits remain isotropic. This is because the BH binary kicked out stars with radial 
orbits more efficiently than those with circular orbits in inner region. 

For run MIT, (3 around the center is close to zero (isotropic) at T = 3, while that in 
outer region is positive. Thus, that stars are radially anisotropic around the outer edge of 
the central cusp, in the triple BH case. In this model, (3 around the center decreases after 
merging at T = 6 like the two BH case. This is because of the slow evolution of the binary 
BH through the interaction with field stars. 



4. DISCUSSION AND CONCLUSION 
4.1. Possibility to Find Triple BH Systems 



In previous sections, we have seen that the density and velocity structure of galaxies 
with triple BHs is different from that with two BHs. In Paper I, we showed that two of three 
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BHs in a galaxy merge within several dynamical times. The standard hierarchical clustering 
scenario of galaxy formation suggests that galaxies frequently merge and form ellipticals. 
Many elliptical galaxies, therefore, might host a binary BH. It is likely that many of them 
had a triple BH system once upon a time. Furthermore, there may be galaxies with a triple 
BH system. The lifetime of a triple BH system in a galaxy is several times the dynamical 
time of the core of the galaxy, which is about 10 8 years (100 M years). If we assume every 
elliptical galaxy have binary BH, and all large ellipticals are formed through mergings of 
elliptical galaxy. Speaking, 1% of large ellipticals could have triple BH systems now. 

One possible way to discriminate between a galaxy with binary or single BH and that 
have (or had) triple BH is the measurement of the cusp radius and density. Our results 
suggest that the central cusp of a galaxy with three BHs is larger by a factor of few than 
that with binary BH. Thus the galaxy with a large cusp and low density in the central region 
might have (or had) three BHs. 

Another possible way is to use the difference of velocity anisotropy. Anisotropy param- 
eter of galaxy with three BHs is positive (radial) in the outer region due to the heating by 
BHs. 

The anisotropy parameter (3 is estimated by the line of sight velocity and their higher 
moment, assuming that the galaxy is spherical. However, with this method it is difficult to 
obtain local a nisotropy param eter in the central region of galaxies, unless the central density 



cusp is steep dGerhard1ll993h. Another method is to fit the model with observational data 



( jCretton et al. 



2000l ; iGebhardt etUboOolliooi l. Figure 10 in bebhardt et all J^OoJ shows 



anisotropy profile of galaxies. The shapes of some galaxies profile in their paper, for example 
all "weak-cusp" galaxies (NGC3608, NGC4291, NGC4649), are similar to of our result in 
Figure These galaxies might contain or have contained triple BHs. 



4.2. Conclusions 

In this paper, we investigated the effect of two- or three-BH systems on the structure of 
galaxies. We found that if the galaxy contains three BHs, (1) multiple three body scattering 
events reduce the density of the central cusp and increase the cusp radius and (2) in outer 
region, orbits of stars are likely to be radial. These difference allow us to discriminate the 
galaxies with two or three BHs. 



We thank Toshiyuki Fukushige, Yusuke Tsukamoto and Keigo Nitadori for stimulating 
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Fig. 1. — Evolution of semi-major axis (top) and eccentricity (bottom) of the most strongly 
bound pair of BHs for all models. Solid and doted curves indicate the results of models MxT 
(triple BH) and MxB (binary BH) respectively. 
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Fig. 2. — Volume (left) and surface (right) density profiles for models MIT (solid) and M1B 
(dashed). Dotted curves show the initial profile. Solid lines indicate the power low p oc r~ 5 . 
Top panel shows the profiles at T = 3, and bottom shows those at T = 6. 
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Fig. 3. — Evolution of Lagrangian radii, radial distance of BHs from the center of the galaxy, 
and the cusp radius, for models MlT(top panel) and MlB(bottom panel). The four solid 
curves show Lagrange radii (0.1%, 0.5%, 1%, 5% of total mass). Thick dashed, dotted and 
dash-dotted curves show the radial distances from the galactic center for BHs. Thick solid 
curves show the cusp radius. 
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Fig. 4. — Evolution of cusp radii (left) and cusp densities (right) for all models. Top and 
bottom panels show the results for triple and two BH models respectively. Solid, dashed, 
and dotted curves represent the results of models Mix, M2x and M3x. 
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Fig. 5. — Spatial density profiles for all runs at T = 6. Solid lines indicate the power law 
p oc r~ a5 . Top panel shows the results of three-BH models, and bottom shows those of 
two-BH models. 
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Fig. 6. — Time evolution of m^e/ normalized by rriBH,tot f° r all models. Top and bottom 
panel show the results of triple- and two-BH models respectively. 
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Fig. 7. — Mass deficit estimated with MM2001's methods and the real value. For MxB 
runs, m def is calculated at T = 1, and for MxT runs at T = 3. 
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Fig. 8. — Velocity dispersion (top) and line-of-sight velocity (bottom) profiles at T = 6 for 
models MIT and M1B. Solid and dashed curves show three- and two-BH cases respectively, 
and dotted curves show the initial profiles. 
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Fig. 9. — Anisotropy parameter for models MIT and M1B. Solid and dashed curves show the 
results of three- and two-BH cases respectively, and dotted curves show the initial profiles. 



